module InletBoundaryConditions

    logical :: InletBCon = .false. ! .true.
    logical :: InletBCreadFile = .false.
    logical :: InitCondInletBCon = .false. ! .true.
    character :: InletBCon_direction = 'r' ! 'x','y','z'

    type tInletBC
        integer :: n
        real(8),allocatable :: RadiusTable(:)
        real(8),allocatable :: UxTable(:)
        real(8),allocatable :: UrTable(:)
        real(8),allocatable :: UtettaTable(:)
        real(8),allocatable :: TempTable(:)
        real(8),allocatable :: PresTable(:)
    end type
    
    type(tInletBC),allocatable :: InletBC(:)

contains

    subroutine CalcNozzleBoundCondInit()
    use Data_cfoam, only: rank
    implicit none
    
        integer :: j, n, OpenFile, ierr
        real(8) :: RadiusTable(16)
        real(8) :: UxTable(16)
        real(8) :: UrTable(16)
        real(8) :: UtettaTable(16)
        real(8) :: TempTable(16)
        real(8) :: PresTable(16)
        real(8) :: Radius, Ux,Ur,Utetta, T, P
        character(100) :: s
        logical :: InletBC_Exists

        ! distribution of boundary inlet parameters 
        data (RadiusTable(j),j=1,16)/ & ! radius jean
                        0, 0.001714543, 0.003851764, 0.007023664,	0.008882004, 0.010368679,&
					    0.012227018,	0.014457029,	0.016315368, 0.018173711, 0.019660382,&
					    0.021518725, 0.023377064, 0.024425507, 0.024955361, 0.025/

    !    data (RadiusTable(j),j= 1,16)/ & ! radius ps90
    !                    0,0.04800719,0.1078494,0.1966626,0.2486961,0.290323,0.3423565,0.4047968,&
    !                    0.4568303,0.5088639,0.5504907,0.6025243,0.6545578,0.6839142,0.6987501,&
    !                    0.7000002/

        data (UxTable(j),j=1,16)/ & ! velocity, x
                        406.7245,420.6993,439.4097,443.9626,448.2266,456.2513,432.3959,333.7299&
                        ,328.3791,332.5684,339.38,349.9492,358.0045,353.8556,299.6601,0/

        data (UrTable(j),j=1,16)/ & ! velocity, r
                        0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0, &
                        0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0/

        data (UtettaTable(j),j=1,16)/ & ! velocity, tetta
                        0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0, &
                        0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0/

        data (TempTable(j),j=1,16)/ & ! temperature
                        559.1461,578.8375,602.9079,600.7621,598.8358,595.0297,501.4955,266.8332,&
                        245.5667,240.5564,237.9458,234.3195,231.6581,233.552,244.145,254.3445/

        data (PresTable(j),j=1,16)/ &
                        31856.41, 31863.25,31734.49,31424.45,31162.23,30900.34,30449.77,29860.01,&
                        29099.75,	28129.51,27170.73,25788.25,	25038.45,25548.08,23896.37,22519.15/

        if(InletBCon.eqv..false.) return
        
        allocate(InletBC(1))
        
        if(InletBCreadFile)then
        
            inquire(FILE='InletBC.dat',exist=InletBC_Exists)
        
            if(InletBC_Exists.eqv..false.)then
                if(rank.eq.0)then
                    write(*,*) '@@ Fatal Error (CalcNozzleBoundCondInit): File InletBC.dat does not exist'
                end if                
                STOP '@@ Fatal Error (CalcNozzleBoundCondInit): File InletBC.dat does not exist'
            end if
        
            open(11,file='InletBC.dat',IOSTAT=OpenFile)
            if (OpenFile.ne.0) then
                if(rank.eq.0)then
	                write(*,*) '@@ Fatal Error (CalcNozzleBoundCondInit): File open error InletBC.dat'
	            end if
	            STOP '@@ Fatal Error (CalcNozzleBoundCondInit): File open error InletBC.dat'
            end if
            
            read(11,*) s
            
            j = 0
            do
                read(11,*,IOSTAT=ierr) Radius, Ux,Ur,Utetta, T, P
                if(ierr.ne.0)then
                    exit
                end if
                j = j + 1
            end do
            
            n = j
            InletBC(1)%n = n

            allocate(InletBC(1)%RadiusTable(n))
            allocate(InletBC(1)%UxTable(n))
            allocate(InletBC(1)%UrTable(n))
            allocate(InletBC(1)%UtettaTable(n))
            allocate(InletBC(1)%TempTable(n))
            allocate(InletBC(1)%PresTable(n))

            rewind(11)
            read(11,*) s
            
            do j = 1,n
                read(11,*,IOSTAT=ierr) InletBC(1)%RadiusTable(j), InletBC(1)%UxTable(j),   InletBC(1)%UrTable(j), &
                                       InletBC(1)%UtettaTable(j), InletBC(1)%TempTable(j), InletBC(1)%PresTable(j)

!                write(1,'(6f25.10)') InletBC(1)%RadiusTable(j), InletBC(1)%UxTable(j),   InletBC(1)%UrTable(j), &
!                                       InletBC(1)%UtettaTable(j), InletBC(1)%TempTable(j), InletBC(1)%PresTable(j)

            end do
            
            close(11)
            
            if(rank.eq.0)then
                write(*,*) 'InletBC.dat has been read successfully!'
            end if
            
!            stop 'hi'
        else
            InletBC(1)%n = 16
            
            allocate(InletBC(1)%RadiusTable(16))
            allocate(InletBC(1)%UxTable(16))
            allocate(InletBC(1)%UrTable(16))
            allocate(InletBC(1)%UtettaTable(16))
            allocate(InletBC(1)%TempTable(16))
            allocate(InletBC(1)%PresTable(16))

            ! distribution of boundary inlet parameters 
            InletBC(1)%RadiusTable(1:16) = RadiusTable(1:16)
            InletBC(1)%UxTable(1:16) = UxTable(1:16)
            InletBC(1)%UrTable(1:16) = UrTable(1:16)
            InletBC(1)%UtettaTable(1:16) = UtettaTable(1:16)
            InletBC(1)%TempTable(1:16) = TempTable(1:16)
            InletBC(1)%PresTable(1:16) = PresTable(1:16)
        end if
        
    
    end subroutine CalcNozzleBoundCondInit

end module InletBoundaryConditions

subroutine CalcNozzleBoundCond(R,U,Temp,Pres)
use InletBoundaryConditions
implicit none
    real(8),intent(in) :: R(3)
    real(8),intent(out) :: U(3),Temp,Pres

    integer :: n,k, i,j
    real(8) :: Radius, Ux,Ur,Utetta
    real(8) :: cosTetta, sinTetta
    ! Tables for parameters for inflow boundary conditions
!    real(8) :: RadiusTable(16)
!    real(8) :: UxTable(16)
!    real(8) :: UrTable(16)
!    real(8) :: UtettaTable(16)
!    real(8) :: TempTable(16)
!    real(8) :: PresTable(16)
    
    n = InletBC(1)%n !16

!    RadiusTable(1:n) = InletBC(1)%RadiusTable(1:n)
!    UxTable(1:n) = InletBC(1)%UxTable(1:n)
!    UrTable(1:n) = InletBC(1)%UrTable(1:n)
!    UtettaTable(1:n) = InletBC(1)%UtettaTable(1:n)
!    TempTable(1:n) = InletBC(1)%TempTable(1:n)
!    PresTable(1:n) = InletBC(1)%PresTable(1:n)


!    ! distribution of boundary inlet parameters 
!    data (RadiusTable(j),j=1,16)/ & ! radius jean
!                    0, 0.001714543, 0.003851764, 0.007023664,	0.008882004, 0.010368679,&
!					0.012227018,	0.014457029,	0.016315368, 0.018173711, 0.019660382,&
!					0.021518725, 0.023377064, 0.024425507, 0.024955361, 0.025/
!
!!    data (RadiusTable(j),j= 1,16)/ & ! radius ps90
!!                    0,0.04800719,0.1078494,0.1966626,0.2486961,0.290323,0.3423565,0.4047968,&
!!                    0.4568303,0.5088639,0.5504907,0.6025243,0.6545578,0.6839142,0.6987501,&
!!                    0.7000002/
!
!    data (UxTable(j),j=1,16)/ & ! velocity, x
!                    406.7245,420.6993,439.4097,443.9626,448.2266,456.2513,432.3959,333.7299&
!                    ,328.3791,332.5684,339.38,349.9492,358.0045,353.8556,299.6601,0/
!
!    data (UrTable(j),j=1,16)/ & ! velocity, r
!                    0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0, &
!                    0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0/
!
!    data (UtettaTable(j),j=1,16)/ & ! velocity, tetta
!                    0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0, &
!                    0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0/
!
!    data (TempTable(j),j=1,16)/ & ! temperature
!                    559.1461,578.8375,602.9079,600.7621,598.8358,595.0297,501.4955,266.8332,&
!                    245.5667,240.5564,237.9458,234.3195,231.6581,233.552,244.145,254.3445/
!
!    data (PresTable(j),j=1,16)/ &
!                    31856.41, 31863.25,31734.49,31424.45,31162.23,30900.34,30449.77,29860.01,&
!                    29099.75,	28129.51,27170.73,25788.25,	25038.45,25548.08,23896.37,22519.15/

! from 0, no negative values
    select case(InletBCon_direction)
    case('r')
        Radius = sqrt(sum(R(2:3)*R(2:3)))
    case('x')
        Radius = R(1)
    case('y')
        Radius = R(2)
    case('z')
        Radius = R(3)
    end select

    i = n   
    do k=1,n
        if (InletBC(1)%RadiusTable(k).ge.Radius)then
            i = k
            exit
        end if
    end do
    
    if(i.eq.1)then
        Radius = 0d0
        i = 2
    end if
    
    if(i.eq.n)then
        Radius = InletBC(1)%RadiusTable(n)
    end if

    Ux = LinInterpolation(InletBC(1)%UxTable)
    Ur = LinInterpolation(InletBC(1)%UrTable)
    Utetta = LinInterpolation(InletBC(1)%UtettaTable)
    Temp = LinInterpolation(InletBC(1)%TempTable)
    Pres = LinInterpolation(InletBC(1)%PresTable)

    U(1) = Ux    
    sinTetta = R(2) / Radius
    cosTetta = R(3) / Radius

    U(2) = Ur * cosTetta - Utetta * sinTetta
    U(3) = Ur * sinTetta + Utetta * cosTetta

contains

    real(8) function LinInterpolation(Table)
    implicit none
        real(8),intent(in) :: Table(n)
        
        LinInterpolation = Table(i-1)+(Table(i)-Table(i-1))/(InletBC(1)%RadiusTable(i)-InletBC(1)%RadiusTable(i-1))* &
                (Radius-InletBC(1)%RadiusTable(i-1))
        return
    end function LinInterpolation

end subroutine CalcNozzleBoundCond
